% Matlab Code for replication of "Optimal Debt Maturity and Firm Investment"
% by J. Jungherr and I. Schott, September 2020
% This code solves the benchmark model with long-term debt & short-term debt 
function solve_model(PARMS)  

% Setup
VAL.period              = 1;
VAL.HISTORY             = cell(1,1);

setup_matrices % Creates empty matrices

% Solve period T-1
solve_T_1

% Iterate forward solution of T-1 in order to generate approximate values for bond price and value function 
simulate_V

% Solve period T-x
for oo = 1:599
    solve_T_x
end

%% Compute model moments

% Compute stationary distribution
stationary_dist

% Compute model moments
model_moments

%% Collect output for Table 2

% 1) Average firm leverage
Tab2.meanlev        = 100*modelmoms.leverage.mean;
% 2) Average credit spread on LTD
Tab2.meanspread     = 100*modelmoms.ltd_spread.mean;
% 3) LTD share
Tab2.LTDshare       = 100*modelmoms.ltd_share.mean;
% 4) Median of firm-level average net investment rate
Tab2.medinv         = 100*modelmoms.i_med_net;
% 5) Median of firm-level s.d. of net investment rate
Tab2.sdinv          = 100*modelmoms.irate_sd_med;
% 6) Total exit rate (exogenous + endogenous)
Tab2.xrate          = 100*moment.xrate;
% 7) Unit mass of firms
Tab2.mass           = moment.mass;

% Save model output
filename            = sprintf('%s%s','Tab2','.mat');
fullname            = fullfile(PARMS.tabfolder, filename);
save(fullname,'Tab2'); % save as .mat file

%% Collect output for Table 3

Tab3.levmean        = 100*modelmoms.leverage.mean;
Tab3.levp25         = 100*modelmoms.leverage.p25;
Tab3.levp50         = 100*modelmoms.leverage.p50;
Tab3.levp75         = 100*modelmoms.leverage.p75;

Tab3.ltdmean        = 100*modelmoms.ltd_share.mean;
Tab3.ltdp25         = 100*modelmoms.ltd_share.p25;
Tab3.ltdp50         = 100*modelmoms.ltd_share.p50;
Tab3.ltdp75         = 100*modelmoms.ltd_share.p75;

Tab3.ltsmean        = 100*modelmoms.ltd_spread.mean;
Tab3.ltsp25         = 100*modelmoms.ltd_spread.p25;
Tab3.ltsp50         = 100*modelmoms.ltd_spread.p50;
Tab3.ltsp75         = 100*modelmoms.ltd_spread.p75;
            
% Save model output
filename            = sprintf('%s%s','Tab3','.mat');
fullname            = fullfile(PARMS.tabfolder, filename);
save(fullname,'Tab3'); % save as .mat file


%% Collect output for Table 5

% 1) Average firm leverage
Tab5_LTD.meanlev        = 100*modelmoms.leverage.mean;
% 2) Average default rate
Tab5_LTD.defrate        = 100*moment.defrate;
% 3) Average credit spread on LTD
Tab5_LTD.meanspread     = 100*modelmoms.ltd_spread.mean;
% 4) Average LTD share
Tab5_LTD.LTDshare       = 100*modelmoms.ltd_share.mean;

% Real variables
Tab5_LTD.AvgK           = moment.AvgK;
Tab5_LTD.AggK           = moment.AggK;
Tab5_LTD.GDP            = moment.GDP;
Tab5_LTD.Cons           = moment.Cons;
    
% Save model output
filename            = sprintf('%s%s','Tab5_LTD','.mat');
fullname            = fullfile(PARMS.tabfolder, filename);
save(fullname,'Tab5_LTD'); % save as .mat file

%% Collect output for Table 6

% 1) Average firm leverage
Tab6_LTD.meanlev        = 100*modelmoms.leverage.mean;
% 2) Average default rate
Tab6_LTD.defrate        = 100*moment.defrate;
% 3) Average issuance-weighted credit spread
Tab6_LTD.meanspread     = 100*moment.iwspread;
% 4) Average LTD share
Tab6_LTD.LTDshare       = 100*modelmoms.ltd_share.mean;

% Real variables
Tab6_LTD.AvgK           = moment.AvgK;
Tab6_LTD.AggK           = moment.AggK;
Tab6_LTD.GDP            = moment.GDP;
Tab6_LTD.Cons           = moment.Cons;
    
% Save model output
filename            = sprintf('%s%s','Tab6_LTD','.mat');
fullname            = fullfile(PARMS.tabfolder, filename);
save(fullname,'Tab6_LTD'); % save as .mat file

%% Collect output for Figure 2

maxz            = 2;

% Compute policy functions over common domain for Figure 2
ltd_pol_plot    = zeros(lenb2,maxz);
cap_pol_plot    = zeros(lenb2,maxz);
lev_pol_plot    = zeros(lenb2,maxz);
ltds_pol_plot   = zeros(lenb2,maxz);
p_pol_plot      = zeros(lenb2,maxz);

% Common debt grid
debtgrid_plot   = linspace(debtgrids(1,1), debtgrids(end,maxz),500);

for zz = 1:maxz
    
    % Find maximum debt level for given Z value
    debt_ss             = find(ltd_pol2(:,zz)*(1-PARMS.GAMMA) < debtgrids(:,zz),1,'first');
    
    ltdF                = griddedInterpolant(PARMS.DEBT{zz},ltd_pol(:,zz),interpmethod,extrapmethod);
    ltd_pol_plot(:,zz)  = ltdF(debtgrid_plot);
    idx                 = ltd_pol_plot(:,zz)*(1-PARMS.GAMMA) > 1.05.*debtgrids(debt_ss,zz);
    ltd_pol_plot(idx,zz) = NaN;
    
    capF                = griddedInterpolant(PARMS.DEBT{zz},cap_pol(:,zz),interpmethod,extrapmethod);
    cap_pol_plot(:,zz)  = capF(debtgrid_plot);
    cap_pol_plot(idx,zz) = NaN;
    
    levF                = griddedInterpolant(PARMS.DEBT{zz},lev_pol(:,zz),interpmethod,extrapmethod);
    lev_pol_plot(:,zz)  = levF(debtgrid_plot);
    lev_pol_plot(idx,zz) = NaN;
    
    % Long-term debt share
    stdF                = griddedInterpolant(PARMS.DEBT{zz},std_pol(:,zz),interpmethod,extrapmethod);
    std_tmp             = stdF(debtgrid_plot);
    totaldebt_plot      = (1+PARMS.C)/(1+PARMS.R)*std_tmp'  + ...
        (PARMS.GAMMA+PARMS.C)/(PARMS.GAMMA+(PARMS.R))*ltd_pol_plot(:,zz);
    ltds_pol_plot(:,zz) = 1./totaldebt_plot.*(PARMS.GAMMA+PARMS.C)....
        /(PARMS.GAMMA+PARMS.R).*((1-PARMS.GAMMA)/(1+PARMS.R))^4.*ltd_pol_plot(:,zz);
    ltds_pol_plot(idx,zz) = NaN;
    
    plF                 = griddedInterpolant(PARMS.DEBT{zz},p_pol(:,zz),interpmethod,extrapmethod);
    p_pol_plot(:,zz)    = plF(debtgrid_plot);
    p_pol_plot(idx,zz) = NaN;
        
end

Av_B            = debtgrids(:)'*probst / sum(probst); % Average stock of existing debt b
Av_BL           = ltd_pol2(:)'*probst / sum(probst);  % Average choice of LTD
Av_K            = cap_pol2(:)'*probst / sum(probst);  % Average choice of capital

debtplot        = debtgrids(debt_ss,zz) * 1.05;

% Histogram
numbins         = 12 * (maxz - 1);
Bin_grid        = zeros(numbins,1);
Bin             = zeros(numbins,2);
% Cover from zero until highest level of debt in Z (maxz)
edges           = linspace(debtgrids(1,1), debtgrids(end,maxz),numbins);
for zz = 1:maxz
    for i = 2:numbins
        % Indicator = 1 if debt within range
        idx             = (debtgrids(:,zz)>=edges(i-1)).*(debtgrids(:,zz)<edges(i));
        % Compute mass of firms in this bin
        Bin(i,zz)       = sum(dist(:,zz).*idx);
    end
end
% Compute median debt level in that bin --> over all firms
for i = 2:numbins
    idx             = (debtgrids(:,1:maxz)>=edges(i-1)).*(debtgrids(:,1:maxz)<edges(i));
    tmp             = debtgrids(:,1:maxz);
    Bin_grid(i)     = median(tmp(idx==1));
end
% First entry cannot be filled
Bin(1,:)        = [];
Bin_grid(1)     = [];

% Normalize by aggregate mass
Bin             = Bin / sum(dist(:));
% Normalize by aggregate debt
Bin_grid        = Bin_grid / Av_B;

% Collect output
Fig2.debtgrid_plot      = debtgrid_plot;
Fig2.debtplot           = debtplot;
Fig2.Av_B               = Av_B;
Fig2.Av_BL              = Av_BL;
Fig2.Av_K               = Av_K;
Fig2.ltd_pol_plot       = ltd_pol_plot;
Fig2.lev_pol_plot       = lev_pol_plot;
Fig2.cap_pol_plot       = cap_pol_plot;
Fig2.p_pol_plot         = p_pol_plot;
Fig2.ltds_pol_plot      = ltds_pol_plot;
Fig2.Bin_grid           = Bin_grid;
Fig2.Bin                = Bin;

% Save output
filename            = sprintf('%s%s','Fig2','.mat');
fullname            = fullfile(PARMS.figfolder, filename);
save(fullname,'Fig2'); % save as .mat file


